Biostat/Biomath M257
This lecture introduces GPU computing in Julia.
GPUs are ubiquitous in modern computers. Following are NVIDIA GPUs on today’s typical computer systems.
| NVIDIA GPUs | H100 PCIe | RTX 6000 | RTX 5000 |
|---|---|---|---|
 |
 |
 |
|
| Computers | servers, cluster | desktop | laptop |
 |
 |
 |
|
| Main usage | scientific computing | daily work, gaming | daily work |
| Memory | 80 GB | 48 GB | 16 GB |
| Memory bandwidth | 2 TB/sec | 960 GB/sec | 576 GB/sec |
| Number of cores | ??? | ??? | ??? |
| Processor clock | ??? GHz | ??? GHz | ??? GHz |
| Peak DP performance | 26 TFLOPS | ??? TFLOPS | ??? TFLOPS |
| Peak SP performance | 51 TFLOPS | 91.1 TFLOPS | 42.6 TFLOPS |
GPUs contain 1000s of processing cores on a single card; several cards can fit in a desktop PC
Each core carries out the same operations in parallel on different input data – single program, multiple data (SPMD) paradigm
Extremely high arithmetic intensity if one can transfer the data onto and results off of the processors quickly
 |
 |
|---|---|
 |
 |
GPU support by Julia is under active development. Check JuliaGPU for currently available packages.
There are multiple paradigms to program GPU in Julia, depending on the specific hardware.
CUDA is an ecosystem exclusively for Nvidia GPUs. There are extensive CUDA libraries for scientific computing: CuBLAS, CuRAND, CuSparse, CuSolve, CuDNN, …
The CUDA.jl package allows defining arrays on Nvidia GPUs and overloads many common operations.
The AMDGPU.jl package allows defining arrays on AMD GPUs and overloads many common operations.
The Metal.jl package allows defining arrays on Apple Silicon and overloads many common operations.
AppleAccelerate.jl wraps the macOS Accelerate framework, which provides high-performance libraries for linear algebra, signal processing, and image processing on Apple Silicon.
The oneAPI.jl package allows defining arrays on Intel GPUs and overloads many common operations.
I’ll illustrate using Metal.jl on my MacBook Pro running MacOS Sequoia 15.4. It has Apple M2 chip with 38 GPU cores.
versioninfo()Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: macOS (arm64-apple-darwin24.0.0)
CPU: 12 × Apple M2 Max
WORD_SIZE: 64
LLVM: libLLVM-16.0.6 (ORCJIT, apple-m2)
Threads: 8 default, 0 interactive, 4 GC (on 8 virtual cores)
Environment:
JULIA_NUM_THREADS = 8
JULIA_EDITOR = codeLoad packages:
using Pkg
Pkg.activate(pwd())
Pkg.instantiate()
Pkg.status() Activating project at `~/Documents/github.com/ucla-biostat-257/2025spring/slides/09-juliagpu`Status `~/Documents/github.com/ucla-biostat-257/2025spring/slides/09-juliagpu/Project.toml`
[13e28ba4] AppleAccelerate v0.4.0
[6e4b80f9] BenchmarkTools v1.6.0
[bdcacae8] LoopVectorization v0.12.172
[dde4c033] Metal v1.5.1
[37e2e46d] LinearAlgebra v1.11.0using Metal, AppleAccelerate
Metal.versioninfo()macOS 15.4.0, Darwin 24.4.0
Toolchain:
- Julia: 1.11.5
- LLVM: 16.0.6
Julia packages:
- Metal.jl: 1.5.1
- GPUArrays: 11.2.2
- GPUCompiler: 1.3.2
- KernelAbstractions: 0.9.34
- ObjectiveC: 3.4.1
- LLVM: 9.2.0
- LLVMDowngrader_jll: 0.6.0+0
1 device:
- Apple M2 Max (4.000 GiB allocated)using Random
Random.seed!(257)
# generate SP data on CPU
x = rand(Float32, 3, 3)
# transfer data form CPU to GPU
xd = MtlArray(x)3×3 MtlMatrix{Float32, Metal.PrivateStorage}:
0.145793 0.939801 0.479926
0.567772 0.577251 0.81655
0.800538 0.38893 0.914135# generate array on GPU directly
# yd = Metal.ones(3, 3)
yd = MtlArray(ones(Float32, 3, 3))3×3 MtlMatrix{Float32, Metal.PrivateStorage}:
1.0 1.0 1.0
1.0 1.0 1.0
1.0 1.0 1.0# collect data from GPU to CPU
x = collect(xd)3×3 Matrix{Float32}:
0.145793 0.939801 0.479926
0.567772 0.577251 0.81655
0.800538 0.38893 0.914135using BenchmarkTools, LinearAlgebra, Random
Random.seed!(257)
n = 2^14
# on CPU
x = rand(Float32, n, n)
y = rand(Float32, n, n)
z = zeros(Float32, n, n)
# on GPU
xd = MtlArray(x)
yd = MtlArray(y)
zd = MtlArray(z);# SP matrix dot product on GPU: tr(X'Y)
# why are there allocations?
bm_gpu = @benchmark Metal.@sync dot($xd, $yd)BenchmarkTools.Trial: 659 samples with 1 evaluation per sample. Range (min … max): 7.475 ms … 8.975 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 7.576 ms ┊ GC (median): 0.00% Time (mean ± σ): 7.591 ms ± 99.875 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▂▄█▆▇▅▅▂▂▂▅▂▃▂▁ ▂▄▂ ▄▆▄▆█████████████████▇████▇▇▆▇▄▅▄▃▄▃▃▄▂▂▃▂▂▄▂▁▁▁▃▁▁▂▁▁▂▁▁▂ ▄ 7.47 ms Histogram: frequency by time 7.85 ms < Memory estimate: 21.24 KiB, allocs estimate: 837.
# SP matrix dot product on CPU: tr(X'Y)
bm_cpu = @benchmark dot($x, $y)BenchmarkTools.Trial: 141 samples with 1 evaluation per sample. Range (min … max): 34.021 ms … 40.975 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 34.774 ms ┊ GC (median): 0.00% Time (mean ± σ): 35.469 ms ± 1.405 ms ┊ GC (mean ± σ): 0.00% ± 0.00% ▅█▂ ▄▅███▇▅▄▅▃▅▁▃▃▅▁▅▃▃▅▄▁▁▃▄▃▄▃▃▅▆▄▅▁▁▁▁▃▁▁▁▁▁▁▁▁▃▁▁▃▁▁▁▁▁▁▁▁▃ ▃ 34 ms Histogram: frequency by time 40.6 ms < Memory estimate: 0 bytes, allocs estimate: 0.
# speedup
median(bm_cpu.times) / median(bm_gpu.times)4.590047354266009# SP broadcast on GPU: z .= x .* y
# why is there allocation?
bm_gpu = @benchmark Metal.@sync $zd .= $xd .* $ydBenchmarkTools.Trial: 547 samples with 1 evaluation per sample. Range (min … max): 8.761 ms … 11.350 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 9.051 ms ┊ GC (median): 0.00% Time (mean ± σ): 9.135 ms ± 347.959 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▄▇█▇█▁▁▂▁▅▃▂ ▃▆████████████▇▆▆▆▇▄▄▄▃▃▃▃▃▁▃▄▃▁▂▃▁▃▂▂▂▁▃▁▃▁▂▃▁▂▃▁▂▁▁▁▂▁▂▁▂ ▄ 8.76 ms Histogram: frequency by time 10.6 ms < Memory estimate: 4.53 KiB, allocs estimate: 177.
# SP broadcast on CPU: z .= x .* y
bm_cpu = @benchmark $z .= $x .* $yBenchmarkTools.Trial: 139 samples with 1 evaluation per sample. Range (min … max): 34.889 ms … 45.114 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 35.275 ms ┊ GC (median): 0.00% Time (mean ± σ): 35.983 ms ± 1.944 ms ┊ GC (mean ± σ): 0.00% ± 0.00% ▇█▇▄▂▁ ▁ ██████▆▆█▁▁▅█▅▁▁▅▁█▁▅▁▅▅▁▁▁▁▁▅▁▁▁▁▁▅▁▁▁▁▅▅▁▅▁▁▁▁▁▅▅▁▁▁▁▁▁▁▅ ▅ 34.9 ms Histogram: log(frequency) by time 44.8 ms < Memory estimate: 0 bytes, allocs estimate: 0.
# speedup
median(bm_cpu.times) / median(bm_gpu.times)3.8975395957005157# SP matrix multiplication on GPU
bm_gpu = @benchmark Metal.@sync mul!($zd, $xd, $yd)BenchmarkTools.Trial: 6 samples with 1 evaluation per sample. Range (min … max): 918.484 ms … 1.044 s ┊ GC (min … max): 0.00% … 0.00% Time (median): 983.040 ms ┊ GC (median): 0.00% Time (mean ± σ): 978.952 ms ± 48.398 ms ┊ GC (mean ± σ): 0.00% ± 0.00% █ █ █ █ █ █ █▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁ 918 ms Histogram: frequency by time 1.04 s < Memory estimate: 960 bytes, allocs estimate: 54.
For this problem size on this machine, we see GPU achieves a staggering 9 TFLOPS throughput with single precision!
# SP throughput on GPU
(2n^3) / (minimum(bm_gpu.times) / 1e9)9.576750193921969e12# SP matrix multiplication on CPU
bm_cpu = @benchmark mul!($z, $x, $y)BenchmarkTools.Trial: 2 samples with 1 evaluation per sample. Range (min … max): 3.960 s … 4.020 s ┊ GC (min … max): 0.00% … 0.00% Time (median): 3.990 s ┊ GC (median): 0.00% Time (mean ± σ): 3.990 s ± 42.588 ms ┊ GC (mean ± σ): 0.00% ± 0.00% █ █ █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁ 3.96 s Histogram: frequency by time 4.02 s < Memory estimate: 0 bytes, allocs estimate: 0.
# SP throughput on CPU
(2n^3) / (minimum(bm_cpu.times) / 1e9)2.221173865839729e12We see >10x speedup by GPUs in this matrix multiplication example.
# cholesky on Gram matrix
# This one doesn't seem to work on Apple M2 chip yet
xtxd = xd'xd + I
@benchmark Metal.@sync cholesky($(xtxd))ArgumentError: ArgumentError: cannot take the CPU address of a MtlMatrix{Float32, Metal.PrivateStorage}
ArgumentError: cannot take the CPU address of a MtlMatrix{Float32, Metal.PrivateStorage}
Stacktrace:
[1] unsafe_convert(::Type{Ptr{Float32}}, x::MtlMatrix{Float32, Metal.PrivateStorage})
@ Metal ~/.julia/packages/Metal/N2ABH/src/array.jl:266
[2] potrf!(uplo::Char, A::MtlMatrix{Float32, Metal.PrivateStorage})
@ LinearAlgebra.LAPACK ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/lapack.jl:3293
[3] _chol!
@ ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:187 [inlined]
[4] #cholesky!#163
@ ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:268 [inlined]
[5] cholesky!
@ ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:267 [inlined]
[6] cholesky!(A::MtlMatrix{Float32, Metal.PrivateStorage}, ::NoPivot; check::Bool)
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:301
[7] cholesky! (repeats 2 times)
@ ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:295 [inlined]
[8] _cholesky
@ ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:411 [inlined]
[9] cholesky (repeats 2 times)
@ ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:401 [inlined]
[10] macro expansion
@ ~/.julia/packages/Metal/N2ABH/src/utilities.jl:10 [inlined]
[11] var"##core#282"(xtxd#281::MtlMatrix{Float32, Metal.PrivateStorage})
@ Main ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:598
[12] var"##sample#283"(::Tuple{MtlMatrix{Float32, Metal.PrivateStorage}}, __params::BenchmarkTools.Parameters)
@ Main ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:607
[13] _lineartrial(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters; maxevals::Int64, kwargs::@Kwargs{})
@ BenchmarkTools ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:186
[14] _lineartrial(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters)
@ BenchmarkTools ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:181
[15] #invokelatest#2
@ ./essentials.jl:1055 [inlined]
[16] invokelatest
@ ./essentials.jl:1052 [inlined]
[17] #lineartrial#46
@ ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:51 [inlined]
[18] lineartrial
@ ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:50 [inlined]
[19] tune!(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters; progressid::Nothing, nleaves::Float64, ndone::Float64, verbose::Bool, pad::String, kwargs::@Kwargs{})
@ BenchmarkTools ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:299
[20] tune!
@ ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:288 [inlined]
[21] tune!(b::BenchmarkTools.Benchmark)
@ BenchmarkTools ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:288
[22] top-level scope
@ ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:461xtx = collect(xtxd)
@benchmark LinearAlgebra.cholesky($(Symmetric(xtx)))BenchmarkTools.Trial: 1 sample with 1 evaluation per sample. Single result which took 7.707 s (0.00% GC) to evaluate, with a memory estimate of 1.00 GiB, over 3 allocations.
@benchmark AppleAccelerate.cholesky($(Symmetric(xtx)))BenchmarkTools.Trial: 1 sample with 1 evaluation per sample. Single result which took 7.550 s (0.01% GC) to evaluate, with a memory estimate of 1.00 GiB, over 3 allocations.
We don’t see GPU speedup of Cholesky at the moment.
# elementwise function on GPU arrays
fill!(yd, 1)
bm_gpu = @benchmark Metal.@sync $zd .= log.($yd .+ sin.($xd))
bm_gpuBenchmarkTools.Trial: 562 samples with 1 evaluation per sample. Range (min … max): 8.789 ms … 9.483 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 8.879 ms ┊ GC (median): 0.00% Time (mean ± σ): 8.896 ms ± 82.151 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▁ ▂▂▂█▃▄▆▅▂▇▄ ▄▂▂▁▃ ▄███████████████████▇▆▆▆▄▅▅▃▄▄▇▅▄▄▂▄▃▂▁▂▂▂▂▂▁▂▁▂▃▁▁▂▁▂▁▁▁▂ ▄ 8.79 ms Histogram: frequency by time 9.2 ms < Memory estimate: 4.53 KiB, allocs estimate: 177.
# elementwise function on CPU arrays
x, y, z = collect(xd), collect(yd), collect(zd)
bm_cpu = @benchmark $z .= log.($y .+ sin.($x))
bm_cpuBenchmarkTools.Trial: 2 samples with 1 evaluation per sample. Range (min … max): 2.721 s … 2.749 s ┊ GC (min … max): 0.00% … 0.00% Time (median): 2.735 s ┊ GC (median): 0.00% Time (mean ± σ): 2.735 s ± 19.800 ms ┊ GC (mean ± σ): 0.00% ± 0.00% █ █ █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁ 2.72 s Histogram: frequency by time 2.75 s < Memory estimate: 0 bytes, allocs estimate: 0.
# Speed up
median(bm_cpu.times) / median(bm_gpu.times)307.9922734440778GPU brings great speedup (>50x) to the massive evaluation of elementary math functions.